#include <UnitTest++.h>
#include "config.h"
#include "grid.h"
#include "optical.h"

using namespace mcrx;
using namespace std;

typedef absorber<double>  T_data;
typedef adaptive_grid<T_data> T_grid;
typedef T_grid::T_grid T_grid_impl;
typedef typename T_grid::T_cell T_cell;
typedef typename T_grid::T_cell_tracker T_cell_tracker;
typedef typename T_grid::T_qpoint T_qpoint;
typedef typename T_cell_tracker::T_code T_code;

// test that column density calculation is correct. this requires a
// "real" grid
struct column_density_fixture1 {
  column_density_fixture1() :
    min(0,0,0), max(1,1,1), n(1,1,1), g(min, max) {

    c=g.begin();

    T_densities rho(1); rho=1;
    c.cell()->set_data(new T_data(rho, vec3d(0,0,0)));

    dir=vec3d(0,1,0);
    dir/=sqrt(dot(dir,dir));
    r.set_position(vec3d(0.5,0.5,0.5));
    r.set_direction(dir);
  }

  vec3d min;
  vec3d max;
  vec3i n;

  T_grid g;
  T_cell_tracker c;

  vec3d dir;
  ray_base r;
};

TEST_FIXTURE(column_density_fixture1, columndensity1)
{
  T_densities n(1);
  bool hit_max;
  T_float len = 
    c.intersection_from_within(r,blitz::huge(T_float()), n, hit_max);

  CHECK_CLOSE(0.5,len,1e-10);
  CHECK_EQUAL(false, hit_max);
  CHECK_CLOSE(0.5,n(0),1e-10);
  CHECK(c.at_end());
}

struct column_density_fixture2 {
  column_density_fixture2() :
    min(0,0,0), max(1,1,1), n(1,1,1), g(min, max) {

    c=g.begin();
    c.cell()->refine();
    c.dfs();
    c.dfs();
    c.cell()->refine();

    c=g.begin();
    const T_cell_tracker e(g.end());
    for(;c!=g.end(); ++c) {
      T_densities rho(1); rho=1;
      c.cell()->set_data(new T_data(rho, vec3d(0,0,0)));
    }

    dir=vec3d(1,.1,0);
    dir/=sqrt(dot(dir,dir));
    r.set_position(vec3d(0,0.25,0.26));
    r.set_direction(dir);
  }

  vec3d min;
  vec3d max;
  vec3i n;

  T_grid g;
  T_cell_tracker c;

  vec3d dir;
  ray_base r;
};

TEST_FIXTURE(column_density_fixture2, columndensity2)
{
  T_densities n(1);
  bool hit_max;
  c=g.locate(r.position(), 0, true);
  CHECK_EQUAL(T_code(0,1), c.code());
  T_float len = 
    c.intersection_from_within(r,blitz::huge(T_float()), n, hit_max);

  CHECK_CLOSE(.5*r.inverse_direction()[0],len,1e-10);
  CHECK_EQUAL(false, hit_max);
  CHECK_CLOSE(.5*r.inverse_direction()[0],n(0),1e-10);
  CHECK_EQUAL(T_qpoint(2,1,1,2), c.qpos());
  r.propagate(len);

  len = 
    c.intersection_from_within(r,blitz::huge(T_float()), n, hit_max);

  CHECK_CLOSE(.25*r.inverse_direction()[0],len,1e-10);
  CHECK_EQUAL(false, hit_max);
  CHECK_CLOSE(.25*r.inverse_direction()[0],n(0),1e-10);
  CHECK_EQUAL(T_qpoint(3,1,1,2), c.qpos());
  r.propagate(len);

  len = 
    c.intersection_from_within(r,blitz::huge(T_float()), n, hit_max);

  CHECK_CLOSE(.25*r.inverse_direction()[0],len,1e-10);
  CHECK_EQUAL(false, hit_max);
  CHECK_CLOSE(.25*r.inverse_direction()[0],n(0),1e-10);
  CHECK(c.at_end());
}

TEST_FIXTURE(column_density_fixture2, max_len)
{
  T_densities n(1);
  bool hit_max;
  c=g.locate(r.position(), 0, true);
  CHECK_EQUAL(T_code(0,1), c.code());
  T_float len = 
    c.intersection_from_within(r,.25, n, hit_max);

  CHECK_CLOSE(.25,len,1e-10);
  CHECK_EQUAL(true, hit_max);
  CHECK_CLOSE(.25,n(0),1e-10);
  CHECK(c.at_end());
  CHECK(!c);
}

struct column_density_fixture3 {
  column_density_fixture3() :
    min(0,0,0), max(1,1,1), n(1,1,1), g(min, max) {

    c=g.begin();
    c.cell()->refine();
    c.dfs();
    c.dfs();
    c.cell()->refine();

    c=g.begin();
    const T_cell_tracker e(g.end());
    T_float rho=1;
    for(;c!=g.end(); ++c) {
      T_densities r(1); r=rho;
      c.cell()->set_data(new T_data(r, vec3d(0,0,0)));
      //cout << c.qpos() << '\t' << rho << endl;
      rho+=1;
    }

    dir=vec3d(1,.1,0);
    dir/=sqrt(dot(dir,dir));
    r.set_position(vec3d(0,0.25,0.26));
    r.set_direction(dir);
  }

  vec3d min;
  vec3d max;
  vec3i n;

  T_grid g;
  T_cell_tracker c;

  vec3d dir;
  ray_base r;
};

TEST_FIXTURE(column_density_fixture3, columndensity3)
{
  T_densities n(1);
  bool hit_max;
  c=g.locate(r.position(), 0, true);
  CHECK_EQUAL(T_code(0,1), c.code());
  T_float len = 
    c.intersection_from_within(r,blitz::huge(T_float()), n, hit_max);

  CHECK_CLOSE(.5*r.inverse_direction()[0],len,1e-10);
  CHECK_EQUAL(false, hit_max);
  CHECK_CLOSE(1*.5*r.inverse_direction()[0],n(0),1e-10);
  CHECK_EQUAL(T_qpoint(2,1,1,2), c.qpos());
  r.propagate(len);

  len = 
    c.intersection_from_within(r,blitz::huge(T_float()), n, hit_max);

  CHECK_CLOSE(.25*r.inverse_direction()[0],len,1e-10);
  CHECK_EQUAL(false, hit_max);
  CHECK_CLOSE(6*.25*r.inverse_direction()[0],n(0),1e-10);
  CHECK_EQUAL(T_qpoint(3,1,1,2), c.qpos());
  r.propagate(len);

  len = 
    c.intersection_from_within(r,blitz::huge(T_float()), n, hit_max);

  CHECK_CLOSE(.25*r.inverse_direction()[0],len,1e-10);
  CHECK_EQUAL(false, hit_max);
  CHECK_CLOSE(7*.25*r.inverse_direction()[0],n(0),1e-10);
  CHECK(c.at_end());
}
